Experimental study of random close packed colloidal particles 
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A collection of spherical particles can be packed tightly together into an amorphous packing known 
as "random close packing" (RCP). This structure is of interest as a model for the arrangement of 
molecules in simple liquids and glasses, as well as the arrangement of particles in sand piles. We use 
confocal microscopy to study the arrangement of colloidal particles in an experimentally realized 
RCP state. We image a large volume containing more than 450,000 particles with a resolution 
of each particle position to better than 0.02 particle diameters. While the arrangement of the 
particles satisfies multiple criteria for being random, we also observe a small fraction (less than 
3%) of tiny crystallites (4 particles or fewer). These regions pack slightly better and are thus 
associated with locally higher densities. The structure factor of our sample at long length scales 
is non-zero, 5(0) = 0.049 ± 0.008, suggesting that there are long wavelength density fluctuations 
in our sample. These may be due to polydispersity or tiny crystallites. Our results suggest that 
experimentally realizable RCP systems may be difi'erent from simulated RCP systems, in particular, 
with the presence of these long wavelength density fluctuations. 

PACS numbers: 82.70.-y, 61.20.-p, 64.70.pv, 64.70.kj 
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I. INTRODUCTION 

Dense packings of hard spheres are an important start- 
ing point for the study of simple liquids, metallic glasses, 
colloids, biological systems, and granular matter 
Of particular interest is the densest possible packing that 
still possesses random structure, "random close packing" 
(RCP), which is important for physics and engineering. 
For example, the viscosity of dense particle suspensions 
diverges when the particles approach the RCP state [l^l . 
In a classic experiment, Bernal and Mason obtained the 
volume fraction of RCP 4>rcp ~ 0.637. They com- 
pressed and shook a rubber balloon which was full of 
metal ball bearings for a sufficiently long time to achieve 
maximum density [l[. Scott and Kilgour also reported 
4>RCP ~ 0.637 by pouring balls into a large vibrating 
container Q. Their results were sensitive to the experi- 
mental method, for example both the frequency and am- 
plitude of vibration. Likewise in computer simulations, 
the value of 4'RCP depends on the protocol. 4>rcp is 
between 0.642 and 0.648 with a rate dependent densifi- 
cation algorithm Q, 0.68 with a Monte Carlo methods Q 
and 0.644 with Lubachevsky-Stillinger packing algorithm 
d, Q. AH of these results are for monodisperse spheres, 
in other words, spheres with identical diameters. 

The variety of results for (/)rcp, in addition to being 
due to the method of preparing the RCP state, perhaps 
also comes from the poor definition of RCP |Tl| . The 
phrase "random close packing" is composed of two terms, 
"random" and "close packing," which are inherently in 
conflict with each other. An ideal random state would 
have no correlation between particles, but the constraint 
that particles cannot overlap already diminishes the ran- 
domness of a physical packing. Furthermore, to get a 
close packing the most efficient method is to pack parti- 
cles into a crystalline array, which is highly non-random 
p^ . For example, a random arrangement of spheres 



can be made denser if it partially includes dense crys- 
talline regions, but then it is less random [l^ [l^. In 
2000 Torquato and coworkers proposed "Maximum ran- 
domly jammed (MRJ)" state as a more tight definition 
of RCP. MRJ states are defined as the least locally or- 
dered structures which are also jammed so that no par- 
ticles can move [ll|. A strictly jammed state should be 
incompressible and unshearable [Tsj . while other defini- 
tions of jammed states can involve external forces [l^ 
or experimental time scales [l7| : the latter can involve 
questions of glassiness. Returning to strictly jammed 
states, one method of quantifying jamming is by con- 
sidering the isothermal compressibility Kt, which is de- 
termined by the structure factor at wave number q = 0, 
Kt ~ ^/p{dp/dp) = S{0)/pkBT where p, p, ks and 
T are density of the material, pressure, Boltzmann con- 
stant, and temperature, respectively. Thus, a strictly 
jammed state requires Kt = 5(0) = since this state 
should be incompressible. Indeed, prior simulation works 
for the strict jammed state of hard spheres show S'(O) « 
to within numerical resolution 0, [H, [l^l • The observa- 
tion S{q — > 0) — > has been termed "hyperuniformity," 
in that the density looks increasingly uniform when con- 
sidered on longer length scales 0. 

The first physics study of the internal structure of a 
random closed packed system that we are aware of is the 
work of Smith, Foote, and Busang [l^. In 1929, they 
studied the packing of shot and used acid to mark the 
contacts between spheres, reporting the contact numbers 
for 1,562 particles taken from the interior of a sample 
with 2,400 particles. In the 1960's, Bernal first studied 
500 particles taken from the interior of a sample with 
5,000 particles m, and later 1,000 particles [13 • In more 
recent times, 16,000 spheres were studied by Slotterback 
et al. who used an index-matching fluid and laser-sheet 
illumination to find the positions [2l| . Aste et al. used x- 
ray tomography to study several different granular pack- 
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ings containing 90,000 particles in an interior region [22j ]. 
These experiments provide useful data for testing theo- 
ries and studying the properties of RCP packings on large 
length scales. 

In this article, we use a sedimentcd dense colloidal sus- 
pension as an experimental realization of a RCP mate- 
rial, in the loose sense of RCP rather than the strict 
sense of a MRJ state. We study the detailed structure 
of our sample with confocal microscopy, which can de- 
termine the three-dimensional positions of the particle 
centers to high accuracy. By carefully imaging overlap- 
ping regions, we observe a large volume containing over 
450,000 particles. Our data set is available online [23| . 
The data are used to determine which features of our re- 
alistic RCP system are similar to the stricter ideal MRJ 
packing. The sample satisfies several criteria for random- 
ness, for example, having only a tiny fraction of particles 
having even short-range order. However, in contrast to 
simulated MRJ packings, we find the isothermal com- 
pressibility is not zero, thus suggesting that in at least 
this particular experimental realization of a RCP system, 
there are differences with simulations. 

It is important to note that our colloidal experiment 
differs in several particulars from both granular experi- 
ments (such as the early ones with ball bearings [ll, Q ) 
and simulations. First, the particles are not all identical; 
they have a polydispcrsity of 5% in their diameters. Sec- 
ond, as the RCP state is formed by sedimentation, the 
particles have a chance to diffuse due to Brownian mo- 
tion. In some situations this motion can help particles re- 
arrange into crystalline packings, if the sample has a vol- 
ume fraction in the range 0.49 < (j) < 0.58 [Ij,!!!!. While 
our experimental preparation method avoids full crystal- 
lization, it is plausible that the sample could be more 
ordered as a result of subtle rearrangements as particles 
sediment toward their final positions. However, conven- 
tionally such sedimentcd colloidal samples are thought 
of as RCP states. Our primary motivation is to use our 
sample to discern properties of the RCP state, and test 
the applicability of ideas derived from simulation. 



II. METHOD 

A. Sample preparation 

We use poly(methyl methacrylate) (PMMA) parti- 
cles sterically stabilized with poly-12-hydroxystearic acid 
[2^. To visualize the particles, they are dyed with rho- 
damine 6G [13]. The mean diameter d of our particles is 
d = 2.53 fim with an uncertainty 1%. Additionally the 
particles have a polydispcrsity of ^ 5%. According to 
prior simulations, the volume fraction for random close 
packing (j)RCP is between 0.64 and 0.66 for a 5% polydis- 
perse system [28l - l30| . which is almost same as (j)RCP for 
monodisperse spheres. References H, [111, H^l point out 
that the specific value often depends on the simulation 
details. 



We use a fast laser scanning confocal microscope (VT- 
Eye, Visitech) which yields clear images deep inside our 
dense samples. Despite the high density, the colloidal 
particles can be easily discerned as shown in Fig. [1] We 
acquire three-dimensional (3D) scans of our sample yield- 
ing a 62.7 X 65.4 X 30 ^m"^ observation volume for each 
image. As the sample is jammed, particles do not move 
and we can scan slowly to achieve very clean images: each 
3D scan takes about 30 s. Within each 3D image, parti- 
cles are identified within 0.03 /im in x and y, and within 
0.05 /xm in z [23,|3l|. 




FIG. 1: Confocal micrograph of the colloidal sediment in (x, 
y) plane. The image was taken 30 /im inside the sample. The 
scale bar represents 10 ^m. The arrow indicates the direction 
of gravity during sedimentation. 

The PMMA particles are initially suspended in a 
mixture of 85% cyclohexylbromide and 15% decalin by 
weight. This mixture closely matches both the density 
and refractive index of the particles [27l |. Then, to in- 
duce the particles to sediment, we add a small amount of 
decalin to slightly decrease the density of the dispersion 
fluid. 

We can quantify the importance of sedimentation by 
computing the nondimensional Peclet number. This is 
the ratio of the time for a particle to diffuse its own ra- 
dius d/2 to the time for it to fall a distance (i/2. The 
diffusion time scale is td = d^/(8Z)), using the diffu- 
sion constant D, which for our particles and solvent is 
Z) = 0.1 fjjr? /s. This gives us r^i = 6 s. The sedimenta- 
tion time scale is rg = d/{2vsed)- We observe the height 
of the sediment as a function of time in a macroscopic 
sample of dilute particles, and find fgod — 0.035 fim/s, 
giving us Ts = 32 s. Thus Pe w 0.2, suggesting that par- 
ticles can diffuse long distances while they sediment; an 
alternate implication is that hydrodynamic interactions 
between particles due to sedimentation are perhaps less 
important than diffusion [s^. Prior to the start of sed- 
imentation, the initial volume fraction is about 0.30 (in 
the stable liquid phase). We stir the particles by ultra- 
sonic wave before sedimentation to avoid the Rayleigh- 
Taylor instability [33| . During sedimentation, the sample 
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passes through the volume fraction range where crystals 
can be nucleated, approximately 0.51 < (j) < 0.60 for 
our sample with 5% polydispersity [13, [SI] . Wc do not 
observe crystals in our final sample, and the most likely 
explanation is that sedimentation happens faster than 
nucleation, which is quite slow for polydisperse samples 
[3^ [sgI [37| . For samples with Pe < 0.1, we do observe 
crystallization, although we have not carefully studied 
the critical Pe for which crystallization is suppressed; 
see Ref. [S^l for further discussion. 

Given that diffusion is faster than sedimentation, the 
sample readily equilibrates, at least at low volume frac- 
tions as the sedimentation starts. Hence, we believe that 
our final state is well-defined and insensitive to the initial 
state. During sedimentation, the Stokes drag force act- 
ing on the particles is given by = 'i-Krjdv, with viscosity 
?7 — 2.18 mPa-s and v = t'sod- The buoyant weight of the 
particles is given by Wh = ApircPg/G with g the acceler- 
ation due to gravity and Ap being the density difference. 
Balancing the gravitational force with the drag force, we 
can estimate the density difference as Ap = 0.038 g/cm^. 
For reference, the particle density is 1.2340 g/cm^. Bal- 
ancing the gravitational energy Wbh with the thermal en- 
ergy ksT lets us determine the scale height h = 1.8 ^m 
(using fcs as Boltzmann's constant). The small scale 
height suggests that in the final sedimented state, there 
will be no density stratification except right at the in- 
terface between the dense sediment and the remaining 
solvent; that interface will have a thickness 0(h). 

During the sedimentation process, it takes about 2 
days for the sample to initially sediment to the bottom 
and form a glassy state. However, the sedimentation 
speed is slow at high (j) [H, [3§|. Thus, we wait 90 days to 
complete the sedimentation before we put the sample on 
the microscope. We also re-checked the sample 300 days 
after the initial sediment, and found the same results as 
a 90 day old sample. 

We use the convention that the y direction is the axis 
corresponding to gravity during the sedimentation pro- 
cess (see Fig. [I}. The sample chamber is made from 
glass slides and coverslips, sealed with UV-curing epoxy 
(Norland), with the sample dimensions being a; = 6 mm, 
y = 20 mm, and z = 0.14 mm. When we measure the 
structure, we lay our sample on the microscope; that is, 
the optical z axis is parallel to gravity and the microscope 
looks into the thinnest dimension of the sample chamber 
(for ease of viewing). In the highly concentrated sample, 
any subsequent gravity-induced particle rearrangements 
arc much slower than our measuring time. In particular, 
we do not observe any particle flow in the sample, and 
the structure docs not change at all during measurement. 
Near the flat coverslip of our sample chamber, particles 
layer against the wall [1^, [llj. To avoid influence of this, 
we take our 3D images at about 1 mm above the y axis 
sample chamber bottom and at about 15 pro. above the 
glass slide along the z axis. Simulations show that wall 
effects decay fairly rapidly 4 diameters = 10 pm) [40| . 
and in our data we see no density fluctuations as a func- 



tion of the distance z away from the coverslip. 

Of course, sedimentation with hydrodynamic interac- 
tions and Brownian motion is not a protocol followed 
in simulations of RCP states. The algorithm developed 
by Lubachevsky and Stillinger considers hard particles 
moving ballistically [3. The particles start very small, 
and continue interacting as they gradually are swelled 
until the system jams. The method of O'Hern and co- 
workers is similar, starting with small particles that grow, 
but their particles are not infinitely hard, nor do they 
have velocities [12] . Rather, the simulation proceeds until 
the particles are maximally swelled but non-overlapping, 
thus giving the final hard-sphere state. Tobochnik and 
Chapin devised a similar algorithm which used Monte 
Carlo moves to eliminate overlaps [H]. These "expand 
and eliminate the overlap" methods are similar to an ear- 
lier method due to Jodrey and Tory which slowly shrank 
spheres, sliding pairs of spheres linearly to minimize their 
overlap, until all spheres had no overlaps [1]. These meth- 
ods all have the strength that the RCP state is generated 
isotropically, in contrast to our experiment where gravity 
sets a direction. (As discussed below, we do not see any- 
thing special about the direction of gravity in our data.) 
Our experimental method does have the feature that our 
spheres never overlap, in contrast to algorithms where 
overlaps arc allowed at intermediate stages [3, [1, [13, [12] , 
although it is not obvious that intermediate stage over- 
laps would cause substantially different results in the final 
state. In some ways our experimental protocol is similar 
to a method by Visscher and Bolsterli from 1972 [43j . 
Their algorithm dropped particles at random positions 
until the particles collided with the floor or a previous 
particle; the falling particle then rolls downhill until it 
reaches a locally stable position. In our experiment, all 
the particles fall simultaneously, and also their Brown- 
ian motion gives them the ability to find better packings 
than the Visscher and Bolsterli algorithm. 



B. Connection of 3D images 

To take a large ensemble, we scan a grid of 3D images 
with a small amount of overlap in x and y. Wc compute 
particle positions from each image and then we connect 
one image to an adjacent overlapping image. Particles 
are considered as superimposed when \rij — rik\ < 0.2 
pm where r*y- is the position of particle i in image j 
and fik is the position of particle / in adjacent image 
k. To achieve this, we apply small displacement shifts 
Ax, Ay and Az to one image, and look for the fraction 
/ of superimposed particles within the overlapped zones. 
Figure [2Ua) shows / in a {Ax, Ay) plane with the reso- 
lution of one pixel accuracy, and we find one spot where 
/ ~ 1 . The secondary ring around the central spot corre- 
sponds to the first peak of the pair correlation function, 
where some coincidences between particle positions are 
expected. While Fig. ^a.) shows / in a two-dimensional 
plane, we calculate / using shifts in the z direction as 
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well. We next apply sub-pixel displacement shifts around 
the spot in Fig. [D^a) to better resolve the peak. Finally, 
we calculate the sum of the squared distances between 
the positions of the overlapped particles within each re- 
gion, and find the global choices of shift values that min- 
imize the overall squared error, to provide the most ac- 
curate shift factors for the overlap. Using the shift fac- 
tors, we then link up the particle positions in adjacent 
sections. The coincident particles are replaced by their 
average position. Figure ^h) shows particle positions at 
5 < z < 5.2 /xm after connecting 4 separate overlapping 
images. The particle positions are well-superimposed in 
the overlapping regions. Our sample chamber contains 
approximately 500,000,000 particles. Using the overlap- 
ping image method, we obtain a large 3D data set with 
size 492 fim x 513 fim x 28 fim, containing more than 
500,000 particles. Due to artifacts when identifying par- 
ticles near the image edges, we clip the data evenly from 
the boundaries and our final data set is F = 492 fim x 
513 ^m X 23.5 /im, containing N = 453, 136 particles. 
This gives us (j)RCP = N{Trd^/6)/V = 0.646 ± 0.020, 
with the error bars due to the 1% uncertainty of the 
mean particle diameter. Our value is in agreement with 
simulations that considered polydispersity [H, [l^ . 

We also examine the average number of particles N 
observed as a function of x, y, and z. To do this as a 
function of x, we count the particles which are located 
between x and a; -f 0.2 yum for a sequence of x values; a 
similar procedure is used for N{y) and N(z). The num- 
ber of particles N{x) as a function of x is fairly flat, as 
are N{y) and N{z), as shown in Fig. HJc). However, 
there are small residual oscillations in x with the stan- 
dard deviation of N{x)/{N) being 0.027 and a period of 
approximately 33 /xmw 13d. This is an artifact of our 
connection algorithm, as we connect the images along 
y direction first, then we connect them along the x di- 
rection. If we change this order, we find N(x) becomes 
flat and N{y) undulates. To evaluate effects of this os- 
cillation, we calculate the structure factor and the pair 
correlation function using both connection orders (x first 
or y first), and find almost identical results. Thus, we 
ignore these oscillations. 



C. Detection of ordered particles 

We use a rotationally invariant local bond order pa- 
rameters c?6 to look for crystalline particles j43 - l46j . The 
idea is to calculate for each particle a complex vector 
Qsmii), whose components m depend on the orientation 
of the neighbors of particle i relative to i. Each of the 13 
components of the vector is given by: 
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where iVt, is the number of nearest neighbor particles for 
particle i, fij is the unit vector pointing from particle 
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FIG. 2: (color online) (a) The image plot of the fraction of 
successfully superimposed particles / in a plane of (Aa::, Aj/). 
The dark central region corresponds to / = 1, meaning that 
all possible particle overlaps are successful, (b) The circles, 
triangles, squares and crosses correspond to particle positions 
obtained from 4 separate 3D images, (c) The local number of 
particles observed A'' normalized by the average, as a function 
of each axis after connecting the images. We add an offset to 
the data of y and z so they can be seen clearly. 



i to its jth neighbor, and Yim is a spherical harmonic 
function. The ggm parameters are the coefficients for 
the spherical harmonics in an expansion of the vector 
directions f^, and thus capture a sense of the structure 
around particle i. The I = 6 harmonics are used as it 
is known that on a local level, hexagonal symmetry is 
often present due to packing constraints [ij, |45|. The 
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neighbors of a particle are defined as those with centers 
separated by less than lAd (which is the location of the 
first minimum of the pair correlation function). These 
13-dimensional complex vectors are then normalized to 
unity using 

Then, we compute c?g as: 



TABLE I: The values of Wi {I =4, 6, 8) for ideal structures 
of fee, icosahedron, hep and bee [4J]. We add the notation of 
(i) for each ideal structure. 















We 


Ws 


fcc(i) 


-0.159316 


-0.013161 


0.058454 


icos(i) 




-0.169754 




licp(i) 


0.134097 


-0.012442 


0.051259 


bcc(i) 


0.159317 


0.013161 


-0.058455 



6 

E ' 

m— — 6 



(3) 



deiiyj) is a normalized quantity correlating the local en- 
vironments of neighboring i and j particles. de{i,j) is 
a scalar and its range is —1 < dg^ijj) < 1; d^^i^j) = 1 
would correspond to two particles who have identical lo- 
cal environments, at least identical in the sense captured 
by the ggm data. Two neighboring particles are termed 
"ordered neighbors" if de{i,j) > 0.5. The number of or- 
dered neighbors is decided for each particle. mea- 
sures the amount of similarity of structure around neigh- 
boring particles. NJ^=0 corresponds to random struc- 
ture around particle i, while a large value of NJ^ means 
that particle i and its neighbor particles have similar sur- 
roundings. Following prior work, particles with > 8 
are classed as crystalline particles, and the other particles 
are liquid- like particles |45| . 

We also compute the VF/ parameter to specify local 
structures: face centered cubic (fee), icosahedral struc- 
ture (icos), hexagonal close packed (hep) and body cen- 
tered cubic (bee) [3]. The Wl parameter is defined as 
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Wl = 



rni ,7712 ,1^3 



I I I 

mi m2 TO3 



(4) 

QlmiQlmiQhn^ (5) 



where () corresponds to the average over neighboring par- 
ticles J, mi + TO2 + = 0, and 



3/2 



(6) 



The coefficients 



/ / / 

nil fTT-s 



are Wigner 3j symbols. Similar to the gg™ parameters 
discussed above, the Wi parameters are able to capture 
a sense of the local ordering with /-fold symmetry, and 
have been used before to help classify local structure; see 
Ref . Q . The values of W^ for ideal structures are listed 
in Table U [31 . These ideal structures are unrealistic for 
experimental data, so we generate 50,000 representations 
of each ordered structure and perturb their positions by 



TABLE IL The ranges of values of Wi {I =4, 6, 8) for struc- 
tures with 5% perturbations from ideal structures We 
add the notation of (p) for the perturbed structures. 





Wi 




We 




Ws 




fcc(p) 


-0.085 ~ 


-0.169 


-0.0109 ~ 


-0.0193 


-0.0180 ~ 


0.0640 


icos(p) 


-0.050 ~ 


0.200 


-0.171 ~ 


-0.162 


-0.090 ~ 


0.090 


hcp(p) 


0.067 ~ 


0.138 


-0.036 ~ 


-0.004 


0.000 ~ 


0.080 


bcc(p) 


0.152 ~ 


0.161 


-0.015 ~ 


0.021 


-0.072 ~ 


0.060 



5 % of the particle diameter, to match the polydisper- 
sity of our experimental particle sizes. This gives us a 
distribution of Wf for each ordered structure (Table HIl . 
Within our experimental data, we calculate W^ (/ = 4, 
6, 8) for each particle. A particle is classed as a ordered 
particle when W4, Wq and W% of a particle are simultane- 
ously within the ranges of one structure shown in Table 
im Otherwise, particles are classed as random particles. 



D. Calculating the structure factor 

We compute the structure factor S{q) via a direct 
Fourier transform of the particle position, S{q) = 
7V~^| X]i=i exp(i(f • fi)\'^ where is the particle position. 
S{q) is the average of S{q) over q = q. 

Our large images have two advantages for calculating 
the structure factor S{q). The first is a high resolution 
with respect to q, as the resolution is given by Sq = 2tt/L 
where L is the image size. Our sample size is 492 fj,m x 
513 /xm X 23.5 /xm and this yields Sq = 0.0128 /xm~^. 
The second advantage of a large image is the reduction 
of boundary effects. Our experimental data set does not 
obey periodic boundary conditions, unlike most simula- 
tions. Thus, we need to use a window function to mini- 
mize the influence of the data cutting off at the bound- 
aries, or we need to periodically replicate the data. Both 
these procedures increase S{q) only near q = 0, but it is 
precisely S{q = 0) that is of interest. Larger images allow 
us to go to smaller q with less problems. We checked the 
Hann window, Hamming window, the Blackman window, 
and also using no window function. We find that S{q) 
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varies for q < 0.55 ^m~^, corresponding to qdjl'K = 0.2. 
That is, our results for q > 0.55 ^m~^ are independent 
of our choice of window functions. In what fohows, we do 
not use a window function, and will focus on the results 
for small q but considering only q > 0.55 fim~^. 



a 0.30 



III. RESULTS 



A. Minimal local ordering 



First, we investigate the randomness of our sample. 
We compute the fraction of ordered neighbors in the sed- 
iment of our colloidal suspension using the parameter 
described in Sec. Ill CI Fig. [3ja) shows the probability of 
finding particles with a given number of ordered neigh- 
bors No- Following prior works, particles with iV* > 8 
are classed as crystalline particles, and the other parti- 
cles arc liquid- like particles [45|. We find the fraction of 
crystalline particles is below 0.03, and that these parti- 
cles are well-dispersed throughout the sample, and shown 
in Fig.[3][b). At most, we see small crystallites composed 
of 3 or 4 crystalline particles which arc nearest neigh- 
bors. Furthermore, the fraction of particles which No is 
below 3 is over 0.8. It means that the coordinate particle 
arrangement of over 80% of particles arc not similar to 
those of nearest neighbor particles. The effects on the 
structure by the spatial distribution of crystalline parti- 
cles will be discussed below. We consider that our system 
is essentially randomly packed as the crystalline particles 
are a quite low fraction and well-dispersed. 

Wc also compute the fraction of specific local ordered 
structures: fee, icosahedron structure (icos), hep andbcc. 
The importance of those structures was emphasized over 
50 years ago by Frank foi' example, the icosahedral 
arrangement has a significantly lower energy than an hep 
or fee cluster for simple Lennard- Jones potentials. To 
specify local ordered structure, we compute Wl (1 = 4, 
6, 8) parameters for each particle Q (see Sec. Ill C|) . We 
find that the fraction of particles that are fee, icosahe- 
dron, hep and bee are 0.0020, 0.0001, 0.0066 and 0.0014, 
respectively. The sum of those fraction is ^ 0.01 and 
this is consistent with the result of No analysis. Again, 
this suggests that the sample is randomly packed. In 
addition, it is interesting that icosahedron is the least 
fraction we observe in our packed hard sphere-like par- 
ticles, whereas icosahedral structure is most stable local 
structure for Lennard- Jones potentials [43|. This is con- 
sistent with many prior observations, and recent simula- 
tions suggest that icosahedral structures are indeed not 
as relevant for random close packed spheres as polytetra- 
hedrons are more favored local structures [i^. We also 
find that the fraction of hep ordering in our sample is 
higher than that of fee. 
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FIG. 3: (color online), (a) The probability of the number 
of ordered neighbors A^o. When a particle has A^^ > 8, it is 
classed as crystalline. The circles corresponds to local volume 
fraction calculated from Voronoi cell volume, averaged over 
all particles with No having the given value, (b) The spatial 
distribution of crystalline particles in 3D image. This image 
is 115 fj,m X 115 fim x 23.5 pm. The other particles are not 
drawn, to better show the crystalline particles. 



B. Voronoi cell volume distribution 

Next, we study the local volume fraction of the sed- 
iment at the particle length scale. We compute the 
Voronoi decomposition which is a unique partitioning of 
space. Each particle is within its own Voronoi polyhe- 
dron, and the Voronoi polyhedron is the region of space 
which is closer to the given particle than any other par- 
ticle [i^. We calculate a volume for each Voronoi cell 
except for those cells located on the boundaries, which 
have incorrectly defined volumes. 

We compute the local volume fraction for each particle 
as 4>i = 7r{d)^ /6Vi where 4>i and Vi arc the local volume 
fraction and Voronoi cell volume for particle i, respec- 
tively. We use the average diameter d since we can not 
detect each particle diameter. The circles in Fig. EKa) 
show the average local volume fraction as a function of 
No- We find that the local volume fraction is almost con- 
stant at iVo < 2, but increases with larger No- This result 
means that few highly ordered particles have a higher lo- 
cal volume fraction than random particles. It is natural 
since ordered phase such as fee crystal is the most packed 
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phase and this tendency is suggested by previous reports 

Next, we compute a distribution of Voronoi cell vol- 
ume. Aste and coworkers proposed a universal function 
of the distribution of Voronoi cell volumes [l^] , and the 
form is described as 

where (V) is the average of the Voronoi cell volumes. It is 
worth noting that the only adjustable parameter in Eq. [7| 
is fc, other than Vmm which is constrained, k is termed 
the "shape parameter" and corresponds the number of 
elementary cells composing the Voronoi cell For in- 
stance, the value of A: is 1 in an fee crystal, while k is close 
to the number of nearest neighbor particles (about 12 or 
13) in random structure [l^]. We choose Vmin = 0.694(i'^, 
which is the smallest Voronoi cell that can be built in a 
packing of monodisperse spheres . Figure |4l^a) shows 
a distribution of the Voronoi cell volumes as a function of 
{V — Vmin)l{^) — Vmin)- The shapc of the distribution 
is asymmetric and not Gaussian, that is, the distribution 
is narrow at small volumes and broad at large volumes. 
We fit the distribution of Voronoi cell volumes with Eq. [7] 
and obtain k = 13.1 (the solid line in Fig. 2]). The tail 
of the distribution is broader than the fitting line, per- 
haps due to the 5% polydispersity of our particles. The 
k value was investigated in experiments using small glass 
beads ('^ 250 fim) in water I22l . acrylic spheres with dif- 
ferent preparation methods [23 and larger glass beads 
3 mm) in oil [2l| . Those similar experiments found 
11 < /fc < 13 [m and close to fc = 14.2 ± 0.6 [13 for 
random sphere packing, k varies with each experiment 
since k slightly depends on the polydispersity. Our exper- 
imentally observed value of k = 13.1 is consistent with 
those prior experiments. This is further evidence that 
the arrangement of our sample is random. We note that 
a universality of k value is proposed of the distribution 
of Voronoi cell volumes for random sphere packing, with 
the evidence coming from experimental with non Brown- 
ian particles [2ll . |22| . Our agreement with the prior work 
suggests that our close-packed sediment is not strongly 
different despite the Brownian motion that the particles 
have during sedimentation. 

Within each Voronoi cell, we now consider the posi- 
tions of particles relative to the Voronoi cell "center of 
mass." We compute a vector Ar.i = ri ~ gi where is 
the position vector for particle i and gi is the position 
vector for the center of mass of the Voronoi cell which 
include particle i. Figure UJ^b) shows the distribution of 
each axis component of Ar^ . We find almost all particles 
are located at the centers of their Voronoi cells within the 
resolution of particle tracking {^^ 0.05 /^m), even along 
the direction of gravity. 




0.0 0.5 1.0 1.5 2.0 



iV-V^in)l{<V>-Vmin) 




Ar (nm) 

FIG. 4: (color online) (a) Distribution of the Voronoi cell 
volumes plotted as a function of (V — Vmin)/((V) — Vmin)- 
The solid line is a fitting line with Eq. [7]using k = 13.1. (b) 
Distribution of the position differences between the particle 
position and the center of mass of the Voronoi cell. Almost 
all particles are located at the center of mass of the Voronoi 
cell. The values on the horizontal axis only go from -0.04 to 
+0.04 /^m, much less than the particle diameter d — 2.53 /im. 



C. Density fluctuations 

We next check whether the sediment is in a "strictly 
jammed state" or not. As we mentioned above, S{0) = 
is required in strict jamming states since strict jamming 
states should be incompressible (equivalently, hyperuni- 
form 0)- To obtain S{0) value, we directly calculate S{q) 
from the particle positions. The inset in Fig. [SJa) shows 
a image plot of the structure factor in a plane of (g^, Qy) 
where and Qy are the x and y components of vector q. 
S{q) is quite isotropic even though y is the direction of 
gravity, again further implying our sediment is randomly 
packed. We average S{q) over q ~ \q\ and obtain S{q), 
shown in Fig. [5{a). Figure Efb) shows S{q) near q~Q 
(circles). S{q) increases near q = because of computa- 
tional artifacts (see Sec. Ill D| ): we find that S{q) is reliable 
over qd/2TT > 0.2, indicated by the vertical dashed line 
in the figure. We fit S{q) with a linear function between 
0.2 < qd/2TT < 0.5 and obtain S'(O) = 0.049 ± 0.008 by 
extrapolation. S{q) is also well fitted by a parabolic func- 
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FIG. 5: (color online) (a) The structure factor S{q) as a func- 
tion of wavenumber q. The inset is the quarter image of the 
structure factor in a plane of {q^, qy). (b) The expanded 
view of S{q) (circles) and Sc{q) (squares) near g = 0. S{q) 
increases near g = because of computational artifacts due 
to the finite size of our data set; the data are reliable for 
qd/2n > 0.2, indicated by the vertical dashed line. The solid 
line is a fitting line with a linear function over the data with 
0.2 < gd/27r < 0.5. We obtain 5(0) = 0.049 by interpolating 
the fitting line to g = 0. On the other hand, the structure 
factor for the crystalline particles Sc{q) decreases with larger 
g and it is well fitted by Ornstein-Zernike function (dashed 
line) over 0.2 < gd/27r < 0.6. 



tion between 0.2 < qd/2TT < 0.5 and S{0) is almost the 
same. Our data are insufficiently strong to determine if 
the linear fit or parabolic fit is more reasonable 0. Our 
uncertainty (0.008) is determined by trying the different 
Fourier transform windowing functions, in combination 
with linear or parabolic fits: all possible combinations 
yield values within the range S{0) = 0.049 ± 0.008, and 
thus we state with confidence that 5(0) ^ 0. Donev, 
Stillinger and Torquato obtained 5(0) = 6.1 x 10"'' by 
numerical simulation with one million monodisperse par- 
ticles and our experimental value of 5(0) is about 
100 times larger than simulation result, a significant dif- 
ference well beyond the uncertainty of our data. This 
implies that our sample is much more compressive than 
the structure found by simulation. 

To support this result, we use a real space function: 
the pair correlation function g{r) shown in Fig. [B] g{r) 



1 2 3 4 5 6 

r/d 

FIG. 6: The pair correlation function g{r) of the sediment. 
The inset is enlarged view at 2 < r/d < 6. The solid line is a 
fitting line with Eq. |8]over r/d > 2. 

at r/d > 2 is well fitted by a exponentially damped os- 
cillatory function [H, [13] described as: 

C 

g{r) ^ - cxp(-r/0 cos[A'o(r - m)] + 1 (8) 

where C, ^ and Kq correspond to an amplitude, a charac- 
teristic length of spatial correlation and the period of the 
oscillations, respectively. From the fitting, we obtain C = 
2.27 ± 0.08, i = 1.50d± 0.03d and Kq = 7.55/d±0.01/d 
Again, we compare with the simulation of a strictly 

framed state which yields ^ = l.SSd and Kq = 7.58/d 
. Though Kq is similar between experiment and simu- 
lation, the length scale ^ from our experiment is shorter 
than that of simulation. The decay of g{r) in an ex- 
periment is related to the broadness of each peak, that 
is, g{r) decays quickly when peaks are slightly broad. 
Broad peaks mean that a distance between two particles 
are distributed. Hence, the rapid decay of g{r) are also 
connected with the fiuctuations in density, supporting 
the result 5(0) 7^ 0. Thus, we conclude that the arrange- 
ment in our experiment is not a strictly jammed state. 
It is important to note that uncertainty in particle po- 
sitions will broaden the first peak of g{r), but this docs 
not strongly affect g{r) for larger r as those uncertainties 
do not accumulate over large distances. That is, the true 
separation between two particles has a distance r^j with 
an uncertainty of ±0.06 /xm, which is somewhat signif- 
icant when rij is small and less important when r.^ is 
large. 

There are several possible explanations for this ob- 
served "softness." One possibility is the polydispersity 
of colloidal size which is a crucial reality for experimen- 
tal situations, and a difference with the simulations to 
which we are comparing our data. When particle sizes are 
slightly different, the minimum distance between two par- 
ticles can be changed from (d) and then the first peak of 
g{r) becomes slightly broad, consistent with Fig. [HI This 
small difference adds up over a long distance and then it 
may induce long wavelength fluctuations. In particular. 
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the local fluctuations in composition (slightly more large 
particles or small particles) are coupled to the number 
fluctuations. A polydispersity of 5% such as we have in 
our experiment results in S{q — ?> 0) = 0.04±0.01 based on 
simulations jssj , consistent with our data. Unfortunately, 
we cannot determine the individual particle sizes as the 
resolution of optical microscopy blurs particle images on 
the same scale as the variability of particle size. Fur- 
thermore, images of neighboring particles overlap, again 
due to the finite optical resolution. This makes determin- 
ing individual particle size problematic, and prevents us 
from disentangling the influence of polydispersity from 
our data. 

A second possibility is that our sample is not at ran- 
dom close packing due to friction effects between the par- 
ticles, which is quite important for granular packings. It 
is known that granular packings are often looser than 
RCP, with volume fractions as low as w 0.58, termed 
"random loose packing" [s^. By vibrating the system, 
the packing fraction can be increased, perhaps coming 
close to RCP [ll [131. In. our experiment, particles move 
by Brownian motion, and this may let them find the RCP 
state. Furthermore, the particles are sterically stabilized 
to prevent them from sticking together. In general, fric- 
tion is not a concept that is usually applied to colloidal 
particles. However, we cannot completely rule out the 
possibility of some possible occasional attractive interac- 
tion between our particles which might result in friction- 
like behavior, resulting in a slightly loose packing. Small 
amounts of static friction gave nonzero 5(0) values in 
simulations [l8| . 

A third possibility is based on the No dependence 
of the local volume fraction (Fig. El^a)), that is, parti- 
cles in more ordered local environments are packed bet- 
ter. The crystalline particles, which have high local vol- 
ume fraction, are distributed throughout the sample (see 
Fig. |31Ib))- To quantify the spatial distribution of the 
crystalline particles, we calculate a crystalline structure 

factor Sc{q) as Sc{q) = N^^ \ J2iLi exp(*9-^j)l^ where 
Wi = 1 when i particle is classed as crystalline, other- 
wise Wi = 0. Sc{q) is the average of Sc{q) over q = q. 
The square symbols in Fig. [SJb) correspond to Sc{q) and 
we find that Sdq) can be fitted with Ornstein-Zernike 
function [dashed line in Fig. ^h)] . This fit gives us that 
the typical length scale between the crystalline particles 
is 12.8d ± 1.8d. Thus, the spatial distribution of crys- 
talline particles can induce density fiuctuations with long 
wavelength and it might be another reason for our obser- 
vation that S{0) 7^ 0. It is worth noting that the small 
but nonzero fraction of the crystallites (less than 3%) is 
crucial in this conjecture. It is possible that these tiny 
crystallites are due to Brownian motion during the sedi- 



mentation. We are unaware of any measurements of tiny 
crystallites in simulations of random close packing, al- 
though one recent study of a binary mixture of spheres 
used the same order parameter that we do and found the 
average number of ordered bonds (our No) was small [i^ . 
IV. CONCLUSIONS 

We use confocal microscopy to study both the local 
and long-range structure of a random close packed col- 
loidal suspension. We find that the fraction of crystalline 
particles is at most 3%, and furthermore that almost no 
regions in the sample have icosohedral order (less than 
0.01%). These observations suggest that the sample is 
randomly packed. This is further supported by the dis- 
tribution of Voronoi volumes, which is well fit by a pre- 
diction based on a model of random packing. 

We also compute the static structure factor S{q) and 
find that S{0) = 0.05, in contrast to simulations which 
found S{0) = 6 X lO"'' 0. S{0) is proportional to the 
isothermal compressibility, implying that the simulated 
states are essentially incompressible (to within numerical 
precision), while our experimental sample is compress- 
ible. This may be due to the presence of tiny crystalline 
regions in our sample, which are associated with slightly 
higher local density (and thus give rise to long wavelength 
density fluctuations) . Alternatively, it may be due to the 
polydispersity of particle sizes in the experiment 5%). 
This softness {S{0) ^ 0) is crucial to how the sample 
would respond to an external force, for example, shear 
stress. The viscosity and elasticity of a sample are ex- 
tremely sensitive to density near jamming point [i^lssj. 
Near the jammed state, small fluctuations in density re- 
sult in large fluctuations of viscosity and elasticity, which 
can lead to shear instability or cracking [s^. This sug- 
gests that real- world RCP materials may possess nontriv- 
ially different properties from idealized simulations. Our 
work points to polydispersity and sample preparation as 
the possible origin of these differences, both of which are 
worthy of further exploration in both simulation and ex- 
periment. 
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